Identifying "communities" within energy landscapes 
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Potential energy landscapes can be represented as a network of minima linked by transition states. 
The community structure of such networks has been obtained for a series of small Lennard- Jones 
clusters. This community structure is compared to the concept of funnels in the potential energy 
landscape. Two existing algorithms have been used to find community structure, one involving 
removing edges with high betweenness, the other involving optimization of the modularity. The 
definition of the modularity has been refined, making it more appropriate for networks such as these 
where multiple edges and self-connections are not included. The optimization algorithm has also 
been improved, using Monte Carlo methods with simulated annealing and basin hopping, both often 
used successfully in other optimization problems. In addition to the small clusters, two examples 
with known heterogeneous landscapes, LJ13 with one labelled atom and LJ38, were studied with 
this approach. The network methods found communities that are comparable to those expected 
from landscape analyses. This is particularly interesting since the network model does not take any 
barrier heights or energies of minima into account. For comparison, the network associated with a 
two-dimensional hexagonal lattice is also studied and is found to have high modularity, thus raising 
some questions about the interpretation of the community structure associated with such partitions. 
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I. INTRODUCTION 

Networks have been the focus of much attention in the 
last few years [1-3]. This is partly due to the growth 
of two networks — the internet [4] and the world wide 
web [5]. A diverse range of networks have also been 
studied, for example interaction networks [6, 7] of pro- 
teins or genes, social networks such as the actor network 
[8, 9] and networks of collaborations amongst scientists 
[10]. Initially, the focus was on the common features that 
many of these networks possessed, such as their small- 
world character [11] or whether the degree distribution 
was scale- free [8]. 

The current emphasis is on more detailed properties 
of such networks, in particular to understand how these 
properties differ between different types of networks and 
how these differences reflect the origins of the networks. 
Properties that have been studied include correlations 
between connected nodes [12], community structure [13] 
and possible spatial embedding [14, 15]. The focus here 
is on community structure, and as highlighted in a re- 
cent summary of 'the ten leading questions for network 
research' there are important open questions in this area 
[16]. For example, what is the best way to quantify such 
community structure and what are its origins? 

The networks that we focus on in this work are asso- 
ciated with potential energy landscapes. The potential 
energy of a system is a function of the coordinates of all 
Nat atoms. Plotting this would give a 37Vat-dimensional 
'hyper-surface' that is known as a potential energy land- 
scape (PEL) [17]. Energy landscapes have been the fo- 
cus of much attention in recent years, as scientists have 
sought to understand the behaviour of a system, such as 
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the ability of a protein to fold to its native state or the 
super- Arrhcnius slowing down of the dynamics in fragile 
liquids, in term of features of the landscape. 

Some of the most interesting points on PELs are sta- 
tionary points where the gradient vanishes. The type of 
stationary points that receive the most focus are the min- 
ima since they represent locally stable structures. It can 
be particularly useful to partition the PEL into basins of 
attraction surrounding each minimum on a PEL, where a 
basin of attraction is defined as the set of points for which 
following the steepest descent pathways from those points 
will lead to the same minimum. Pioneering work by Still- 
inger [18] has shown that systems generally spend most 
of their time vibrating in a minimum's basin of attrac- 
tion and occasionally hop to another basin, although at 
sufficiently high temperature this time scale separation 
begins to break down [19]. First-order saddle points (of- 
ten called transition states) are particularly important to 
describe this hopping, since the trajectory of the system 
passes along a transition state valley connecting adjacent 
basins during this process. Thus, the minima, their en- 
ergies, their connectivities and the heights of barriers be- 
tween them control the dynamics of the system, i.e. how 
it moves around the PEL, and hence how the structure 
changes. 

The number of minima increases exponentially with 
the number of atoms in a system. Therefore, for all but 
quite small systems, it will be impossible to completely 
characterize the PEL. Instead, the aim is to obtain a 
statistically accurate representation of the PEL from an 
incomplete sample of minima and transition states. This 
problem has essentially been solved for minima, allowing 
the distribution of minima as a function of energy to be 
obtained [20-22]. Hence, a landscape description of the 
thermodynamics is possible. 

However, the situation is much more complex for the 
dynamics, because not only is the distribution of transi- 
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tion states required, but also their pattern of connectiv- 
ities. Particularly little is known about the latter, and 
the aim of our research programme is to get insights into 
the fundamental organization of the PEL by studying 
the PEL connectivity for small systems for which a com- 
plete characterization can be obtained. To achieve this 
we have applied the tools developed to study complex 
networks to the network of minima linked by transition 
states. It should be noted that this representation only 
describes the connectivity of the PEL, its topology. It 
does not take into account the energies of the minima or 
the heights of the barriers between them, its topography. 

Recently we showed for a series of small Lennard- Jones 
(LJ) clusters that the networks show both a small-world 
behaviour and have a power-law tail to the degree distri- 
bution [23, 24]. In these networks the low-energy min- 
ima with large basins of attraction act as hubs in the 
network. Furthermore, a detailed characterization of a 
variety of topological properties, such as the local clus- 
tering coefficients, degree-degree correlations and mea- 
sures of betweenness showed behaviour similar to other 
scale-free networks. In this paper we extend this work by 
applying algorithms designed to identify communities of 
strongly connected nodes to our PEL networks. Our aim 
is to compare the divisions of the PEL achieved by these 
methods, which are based on purely topological informa- 
tion, to methods from the field of energy landscapes that 
also use topographical information. 



II. LANDSCAPES, FUNNELS AND 
COMMUNITY STRUCTURE 

In social networks, people are often grouped into com- 
munities where there are particularly strong interactions 
between the people in a community, and the people share 
some common characteristics or purpose. Attempts to 
quantify these concepts in network terms usually focus 
on the greater density of edges within communities than 
between them. For example, the; algorithm of Girvan 
and Newman [13] discussed in the following section was 
applied to a network where nodes were people in an or- 
ganization and edges linked those that communicated by 
email [25]; the community structure identified by the al- 
gorithm corresponded to informal social groups within 
the organization, as confirmed by the people involved. 

In the PEL networks studied here, communities would 
indicate groups of minima interconnected by many tran- 
sition states. It would be reasonable to expect conver- 
sions between those minima to occur rapidly whereas 
those between different communities/groups of minima 
would occur more slowly. Community structure would 
then have a strong effect on the dynamics of the system. 
If there is strong community structure then groups of 
minima are segregated from each other, probably making 
it more difficult to sample all of the PEL. However, the 
topography of the PEL is also expected to have a strong 
effect on this. For example, a group of minima might 



have many transition states between them and hence be 
classified as a community in this approach, but if the bar- 
riers are high in energy then they are unlikely to be used 
by the system. 

The relationship between dynamics and community 
structure has been put to use previously in chemistry via 
the master equation approach [17]. The master equation 
describes the rate of flow of the occupation probability 
between minima in terms of a matrix of transition rates 
between each individual minima, calculated using uni- 
molecular rate theory. The eigenvectors of this matrix 
are associated with different dynamical processes. If the 
network has two different communities, the slow mode 
associated with transitions between those communities 
has an eigenvector where the sign of the components de- 
pends upon in which community the minimum is. The 
main difference between that approach and this work is 
that the only information used here is the connectivity. 
In network terms the focus is an unweighted, undirected 
adjacency matrix, rather than a transition matrix where 
the strength of the connections depends on properties 
such as the height of the barrier between the two min- 
ima. Our approach is equivalent to studying the system 
at infinite temperature, where all barriers can be equally 
easily crossed. Interestingly, spectral methods have also 
been used to divide networks into communities based on 
the eigenvalues and eigenvectors of the adjacency matrix 
[26, 27]. 

Simple PELs such as those of the small clusters stiidied 
here are shaped like funnels leading down to the global 
minimum [28, 29] and therefore may be expected to have 
no community structure. Other systems, e.g. the 38-atom 
LJ cluster, have multiple funnel structures with large bar- 
riers between the different funnels [29] . In L J38 the global 
minimum is a truncated octahedron and is at the bottom 
of a narrow funnel of face-centred-cubic (fee) structures, 
whereas most low-energy structures are icosahcdral and 
at the bottom of a second wider funnel. Evidence for 
these two funnels comes from analyses of both the to- 
pography of the PEL and the cluster's dynamics [30]. 
The LJ38 disconnectivity graph, which provides a picto- 
rial representation of the barriers between the minima, 
shows a clear division of the low-energy minima into two 
sets, termed funnels because the energy decreases as one 
steps down towards the low-energy minima at the funnel 
bottom. There are low barriers between minima in the 
same funnel, but a much larger barrier between the two 
funnels. Unsurprisingly given that the energy barrier is a 
major determinant of transition rates, there is a separa- 
tion of time scales between transitions between minima 
within a funnel and transitions between the two funnels. 
There is also a thermodynamic solid-solid transition from 
the fee to icosahedral structures just below the melting 
point [31]. There are more icosahedral structures than 
fee structures, so icosahedral structures are favoured en- 
tropically at higher temperatures, whereas the low energy 
global minimum is favoured at low temperatures. 

LJ38 will be studied here to determine whether network 
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methods based purely on connectivity can find groups 
corresponding to the two funnels. The sample of minima 
and transition states we analysed was previously used to 
study the multiple funnel topography and the thermody- 
namics and dynamics of the solid-solid transitions [30]. 
Unlike the small clusters however, this sample is far from 
complete. The method previously used to identify the 
minima in the two funnels is based on a master equation 
approach. 

Any structure of a cluster will have many permuta- 
tional isomers, structures which are identical except that 
the positions of some atoms have been exchanged. For 
a monoatomic cluster, the number of possible permuta- 
tions of each geometrically distinct minimum is A^at- As 
some of these are related by symmetry this number is re- 
duced to 2Na,t^./o, where o is the order of the point group 
of that minimum. 

In our analyses of the PEL connectivity of L J clusters 
so far we have not made any attempt to distinguish such 
permutational isomers of the same minimum. Instead, 
we have just looked at the connections between geomet- 
rically distinct minima. The PELs of small LJ clusters 
are usually described as having a single-funnel topogra- 
phy. However, if one starts to consider different permu- 
tational isomers the situation can become more complex. 
For example, if for LJ13 one considers one atom to have 
a different mass, there are two physically distinguishable 
isomers of the icosahedral global minimum, one with the 
heavy atom in the centre of the icosahedron and one with 
this atom on the surface. As there is a significant barrier 
for exchange of an atom between the centre and surface of 
the icosahedron, the two isomers of the global minimum 
lie at the bottom of separate funnels [32] . 

III. METHODS 

A. Community structure algorithms 

Identifying communities within a network is a non- 
trivial task, because the number of possible divisions of 
the network is very large, especially since, in general, the 
sizes and number of communities are unknown. A num- 
ber of different algorithms have been proposed to attempt 
this task [26, 27, 33-38]. A recent method proposed by 
Girvan and Newman (the betweenness algorithm) [13, 39] 
relies solely on the connectivity of the nodes and has 
proven successful in various networks. The betweenness 
of an edge is based on the number of shortest paths that 
pass along it. If an edge lies between two communities, 
it is expected to have higher betweenness because lots of 
shortest paths will pass through it that connect nodes in 
the different communities (from the definition of a com- 
munity there will be few inter-community edges for the 
shortest paths to choose from). If the edges between 
communities are removed then the network will be bro- 
ken down into components corresponding to the different 
communities of the original network. After each edge is 



removed the topology of the network changes and the 
betweenness must be recalculated. The algorithm is hi- 
erarchical, edges are removed until the network is broken 
down from one community of nodes into N commu- 
nities of one node each with various community splits in 
between. This can be represented in a dendrogram show- 
ing the various splits in a hierarchical manner, a vertical 
line at any point gives the community split at that point. 

To determine which of the many splits found by this 
algorithm is the best, a quantitative measure of the com- 
munity structure is needed. The fraction of edges lying 
within communities should be large for a good commu- 
nity split. This is denoted by Cc where Cc is the frac- 
tion of edges lying within community c. However, the 
maximal value of J2c found if the whole network is 
considered as one community, in which case all edges lie 
within that community and X]c^c=l- To rectify this, 
Newman and Girvan introduced the modularity Q. To 
calculate Q, the fraction of edges that would lie within 
the same communities in a random network where the 
nodes have the same degree is subtracted from J2c^c- 
Therefore, if there are more edges in communities than 
would be expected due to the degree distribution, Q is 
large. A plot of Q against the number of edges removed 
will have a peak with typical values of Q around 0.3-0.7 
if there is community structure. 

The fraction of edges within communities in the real 
network is simply counted, and the expected value for 
a random network is calculated from the degree of the 
nodes. If an edge is picked at random, the probability 
flc that one end of it leads to community c is propor- 
tional to the number of edges that lead there, which in 
turn is simply the sum of the degrees of all nodes in that 
community. Therefore 



The probability that both ends of the edge lead to com- 
munity c is then a^, and so Q is given by 

Q = ^ (e, - al) . (2) 

c 

Having such a measure of community structure also 
naturally leads to a new approach to finding communi- 
ties. Rather than just using Q to evaluate the output of 
a community structure algorithm one can instead try to 
optimize Q directly. In the first application of this ap- 
proach, Newman used a simple agglomerative and greedy 
algorithm [40]. Starting with N communities each con- 
taining one node, communities are joined together one 
node at a time until there is a single community. At 
each step, the two communities that lead to the largest 
increase in Q (or if no increases in Q are possible, the 
smallest decrease) are grouped together. This greedy al- 
gorithm is faster than that based on betweenness, scaling 
as 0{{M + N)N) rather than 0{M'^N). 

Finding the maximal value of Q is a global optimiza- 
tion problem, and while the above greedy algorithm can 
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find community splits with high modularity, there is of 
course no guarantee that it will find the split with max- 
imum modularity. Indeed, given the large search spaces 
for the larger networks we analyse, there is good reason 
to think that it will not find this optimal split. If better 
community splits are sought, more sophisticated global 
optimization algorithms should bo used. 

In this paper a second approach to optimize Q is also 
investigated. This is based on the Monte Carlo method 
with simulated annealing. At each step, a node and a 
community are chosen at random. The community could 
be any of the existing communities, including the one 
that the node is already in, or it could be a new com- 
munity that does not contain any nodes. Moving the 
node from its initial community to the new community 
would change Q by AQ. If AQ is greater than zero then 
the move is accepted, otherwise the move is accepted 
with probability exp(/3A(3). This is the Metropolis cri- 
terion, where (3 represents the inverse temperature. At 
high temperatures many moves are accepted and lots of 
different community splits are sampled, whilst at lower 
temperature fewer splits are sampled but those that are 
generally have higher modularities. The principle of sim- 
ulated annealing involves starting the algorithm at high 
temperature and decreasing the temperature until Q be- 
comes constant because no further moves arc accepted. 
The hope of such an annealing scheme is that the fi- 
nal partition will be the global optimum, but of course it 
could only be a local optimum. To increase the likelihood 
of success, quenches can also be applied periodically. In 
these quenches, AQ is calculated for moving all nodes to 
all communities, and the move with the highest AQ is 
taken. This is repeated until the highest AQ is less than 
or equal to zero implying Q cannot be increased further. 
The simulated annealing method runs quickly, taking a 
similar amount of time as the greedy optimization but 
finding higher values of Q. 

However, simulated annealing is often not a particu- 
larly efficient global optimization method. Indeed, we 
were able to obtain still higher values of the modularity 
using a basin hopping approach that is known to perform 
well in other optimization problems [41]. Basin- hopping 
also uses Monte Carlo, but there are a number of sig- 
nificant differences from the simulated annealing scheme 
described above. Firstly, each step consists of randomly 
changing the communities of a series of nodes, not just 
a single one. Secondly, after each step the new partition 
is quenched and the Metropolis acceptance criterion is 
applied to the values of the modularity for the partitions 
that result from quenching. Then, if a step is accepted, 
the ciirrent partition is updated to that of the quenched 
partition. This algorithm is slower to run to completion, 
but finds high values of Q quickly. 

The initial point for both these approaches can be any 
partition of the nodes into communities. This could be A'' 
communities of one node, but the algorithms are faster 
if the initial partition is either that obtained from the 
greedy algorithm or a community split obtained in a sim- 
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FIG. 1: Variation in the clustering coefficient of the LJ14 
network as edges were rewired to form a random ensemble. 
Data collection began at 5 switches per edge. 

ilar agglomerative manner but using Monte Carlo, which 
is much faster. To create this initial partition, rather 
than checking all possible edges as in the greedy algo- 
rithm, a single edge is chosen at random and added if 
AQ meets the Metropolis criterion. Simulated annealing 
runs were started from the partition obtained from the 
greedy algorithm, which gave much higher values of Q. 
If the algorithm is started from the second more random 
initial partition then the results depend quite strongly 
on parameters such as the cooling rate and frequency of 
quenches. However, for the basin hopping algorithm the 
initial conditions, in general, have only a weak effect on 
Qmax, so this was started from the more random parti- 
tion. This makes basin hopping more feasible for larger 
networks, because the greedy algorithm, although gen- 
erally fast, especially using the algorithm of Ref. [42], 
requires large arrays for large networks. 



B. Randomized networks 

In order to interpret the results of the above algo- 
rithms, it is important to compare the values of Q ob- 
tained to those for an appropriate null model. Given the 
importance of the degree distribution for other network 
properties, the usual null model is an ensemble of ran- 
dom networks that have the same degree distribution as 
the network being analysed. Currently, the most efficient 
algorithm to generate such an ensemble is the rewiring 
method [43, 44]. The starting point is the real network. 
Two edges A-B and C-D are then picked at random and 
rewired to give either A-C and B-D or A-D and B-C. 
Moves arc simply not allowed if they would create mul- 
tiple edges or self-connections. 

As edges were rewired the clustering coefficient was 
measured to follow the randomization. This is shown 
plotted against the number of rewired edges (switches) 
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per edge of the network for LJ14 in Figure 1. All of these 
LJ networks appear to be fully randomized after around 
one switch per edge, so data collection for the ensemble 
began at five switches per edge. It is likely that the clus- 
tering coefficient reaching its equilibrium value indicates 
that the networks are random after one switch per edge, 
but since the algorithm is fairly fast, rewiring more edges 
has no disadvantages. For all random networks with be- 
tween five and ten switches per edge, the probability that 
any two nodes were connected was recorded. 25 random 
networks were also saved at intervals of 0.2 switches per 
edge (roughly the autocorrelation time of the clustering 
coefficient) from the same range. The community struc- 
ture algorithms were then run on this ensemble, recording 
the mean and standard deviation of the modularity. 



IV. RESULTS 
A. LJs-M 

The results of applying the edge betweenness algorithm 
to identify community structure are illustrated in Figure 
2 by the dendrogram and the associated modularity Q 
for LJio (all small cluster networks studied showed qual- 
itatively similar behaviour). 

The dendrogram shows that nodes were removed from 
the initial 'community' of nodes roughly one by one. 
At any stage there is one large community and several 
smaller 'communities' consisting of one or two nodes 
each. The global minimum is the last node remaining 
in the original community, which is perhaps unsurprising 
since it has the highest degree. A plot of the degree of the 
node removed at each stage against the order in which 
they are removed (Figure 3) shows that lowest degree 
nodes are removed first. 

This behaviour has been seen in previous work for ran- 
dom networks without community structure [45] and our 
random networks also have qualitatively similar dendro- 
grams. In the case of a node with degree one, all shortest 
paths to that node must pass along the one edge connect- 
ing it to the rest of the network, making the betweenness 
of that edge 2(A^ — f ) (the factor of 2 appears because in 
calculating the betweenness each path is counted twice) . 
This is true in any network, but the fact that these edges 
have the highest betweenness reflects the relatively low 
betweenness of the rest of the edges in the LJ networks. 
This has been used previously to identify a lack of com- 
munity structure in a network [25, 46]. In that case, when 
enough edges had been removed to make the edges con- 
nected to degree one nodes have the highest betweenness, 
the algorithm was considered to have run to completion. 
Since the magnitude of Q is so small (an order of mag- 
nitude smaller than most networks showing community 
structure) the results seem to indicate that the LJ net- 
works have little community structure. 

Using the greedy algorithm to optimize Q, however, 
does find some community structure, as illustrated by the 
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FIG. 2: (a) Dendrogram showing the results of the between- 
ness algorithm applied to LJio. Each point on the y-axis 
represents a single node in the network (in arbitrary order). 

Drawing a vertical line at any point, as has been done at the 
best splits as determined by Qmax and Q'max , gives the num- 
ber and sizes of communities at that point. In the case of 
the betweenness algorithm, the initial point is the right-hand 
side where all nodes are in a single community. The scale of 
the X-axis is the number of edges remaining in the network, 
(b) Variation in Q (solid line) and Q' (dashed line) as the 
algorithm progresses. 

dendrogram for LJio in Figure 4. The community struc- 
ture is weak compared to other networks — the magnitude 
of Q for the best split, Qmax, is around 0.2-0.3. Qmax 
for the greedy algorithm is much higher than that for the 
betweenness algorithm, therefore the community splits 
from this algorithm are much better. The betweenness 
algorithm has been found to be unsuccessful for dense 
networks (e.g. a food web [13]), since if there are many 
edges between communities the shortest paths between 
the communities can pass through any of these edges and 
the betweenness of any one is unlikely to be high. It is 
noteworthy that in this respect the average degree of the 
LJ networks are significantly larger ((fc) up to 30 [24]) 
than those of many of the networks previously studied. 
The betweenness algorithm has also been able to detect 
communities in weighted networks when it has failed for 
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FIG. 3: Degree of nodes separated from the main community 
by the betweenness algorithm for LJio. 

the equivalent unweighted network [27]. Here we con- 
centrate on improving the quantification of community 
structure and hence the community structure found by 
optimizing this value. 

Q is negative for community splits with many small 
communities (to the left in the dendrogram), implying 
that the community structure is weaker than in a ran- 
dom network with the same degree distribution. The 
reason for this is that in 'communities' of one node there 
can be no edges within that community. However, the 
expression for Q (Equation (2)) uses a predicted fraction 

of edges of = {^i^^ki/2M^ that is clearly greater 
than zero. This effect makes Q smaller when there are 
lots of small communities, as is the case for the between- 
ness algorithm, because the second term in Q predicts 
more edges than are physically possible. The total num- 
ber of edges predicted is equal to the total number of 
edges in the network, so if too many are predicted in one 
area too few must be predicted elsewhere. Therefore, as 
the communities grow this effect is reduced until (5 = 
when there is one large community. Using to pre- 
dict the number of edges in a community if the network 
were random is equivalent to predicting kikj/2M edges 
between two nodes i and j. This number can clearly 
be greater than one for two high degree nodes. For the 
largest network studied, LJ14, the global minimum has 
degree 3201 and there are 61085 edges, so more than one 
edge is predicted between the global minimum and any 
node with degree greater than ~40. This is similar to 
the argument used by Maslov and Sncppcn to explain 
degree-degree anticorrelations seen for the internet [43]. 
The lack of self-connections means Q disfavours small 
communities and the lack of multiple edges means Q dis- 
favours the grouping together of high degree nodes. 

A solution to this problem is to improve the predicted 
fraction of edges within communities. The method used 
here was to create an ensemble of random networks with 
the same degree distribution as the original and with the 
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FIG. 4: (a) Dendrogram showing the results of the greedy 
algorithm applied to LJiq. The initial point for the greedy 
algorithm is at the left-hand side with A'' (64) communities. 
The X-axis represents the number of steps, where each step 
involves joining together two communities. The dcished line 
represents the best split found (with maximum Q). (b) Vari- 
ation of Q as the algorithm progresses. 

constraint that multiple edges and self-connections are 
forbidden. It is then straightforward to calculate the 
probability that there is an edge between any pair of 
nodes in a network taken from that ensemble. The pre- 
dicted fraction of edges within each community in a ran- 
dom network (with the same degree distribution and no 
multiple edges or self-connections) is then the sum of this 
probability over each pair of nodes within the commu- 
nity. This sum (denoted by fc) replaces in Equation 
(2) giving 

c 

It may also be possible to use the analytical method of 
Park and Newman [47] to predict the probability of an 
edge between two nodes. 

The probability of an edge between two nodes is shown 
in Figure 5 plotted against the product of the degree 
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FIG. 5: Average probability of an edge between two nodes 
as a function of the product of the degrees of those nodes, 
calculated for LJ14. Values were calculated from an ensemble 

of random networks and arc used to calculate Q' . The number 
of edges predicted by this method is never larger than 1. The 
straight line shows the value kikj/2M used in the original 
calculation of Q. For the two highest degree nodes, this takes 
the value 33.46. 

of the two nodes for both the original prediction = 
kikj/2M and fij from rewiring, as used in calculating 
Q'. fij is roughly proportional to the product of the de- 
grees of the nodes at cither end of each edge for inecliuin 
degrees, but there is a plateau at higher degrees in fij 
so that fij < 1 for all pairs of nodes, fij is then neces- 
sarily larger than for nodes with lower degree since 
the sum over all pairs of nodes is M. This approach has 
also proved useful for the assortativity coefHcient [12, 24], 
which also requires the predicted number of edges be- 
tween two nodes in a random network. 

Using Q' rather than Q will not affect the between- 
ness algorithm (although it could affect which community 
split is identified as the best), but it will almost certainly 
affect the greedy algorithm. Q' for the betweenness algo- 
rithm is shown in Figure 2 and the dendrogram and Q' 
for the greedy algorithm in Figure 6. 

The overall effect is unchanged; the greedy algorithm 
finds some community structure with a larger value of Q' 
than for the betweenness algorithm, which finds essen- 
tially no community structure. However, the community 
split found by the greedy algorithm is different for the 
two measures of modularity. For the betweenness algo- 
rithm Q' can be compared to Q since both are measured 
for the same community splits. The difference is signif- 
icant for this network, with the ordering of the peaks 
being altered. The highest peak of Q' is further to the 
left because the original Q disfavours small communities. 

Higher values of Q' were found by simulated anneal- 
ing and basin hopping runs. These values are shown 
in Figure 7, along with those for the greedy algorithm 
and the betweenness algorithm, for the real and ran- 
dom networks. The modularity increases roughly log- 
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FIG. 6: (a) Dendrogram showing the results of the greedy al- 
gorithm using Q' applied to LJio. The dashed line represents 
the best split found (with maximum Q'). (b) Variation of Q' 
as the algorithm progresses. 



arithmically with the size of the networks. The rate of 
increase for Q'^^x from the greedy algorithm slows down 
for larger networks, as seen previously [48], simply be- 
cause as the networks become larger the number of pos- 
sible community splits increases and so the probability 
of finding the global maximum of Q' decreases. As the 
networks become larger, they also become more sparse 
(p = 2M/{N{N-1)) decreases). Therefore, there should 
be community splits with less edges between communi- 
ties, and hence greater Q' . 

Q' is greater than zero for the random networks mean- 
ing that it is always possible to find some community 
structure. This feature, coupled with the fact that Q' 
depends on the size and average degree of the network 
[48] , implies that it is important to compare results with 
those for random networks [49]. Q'^^ax several stan- 
dard deviations above the mean of the random networks, 
so although Q^ax fairly low compared to other net- 
works, there is significantly stronger community struc- 
ture than expected if nodes were connected randomly. 
This is perhaps due to the spatial nature of the net- 
works and their high clustering. For example, previous 
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FIG. 7: Q' against N for cluster networks and random en- 
sembles, error bars are one standard deviation from an en- 
semble of 25 random networks. The results are (from bottom 

to top) for the bctwccnncss algorithm, the greedy algorithm, 
simulated annealing and basin hopping. The betweenness al- 
gorithm was only run on clusters with up to 12 atoms (515 
nodes) as it is slow for larger networks. 



work on the contact networks associated with models of 
space-filling disks (specifically the two-dimensional Apol- 
lonian packing) [50] found community structure with 

Qmax = 0.5938. 

Some correlations have been seen between community 
structure, high clustering and assortativity. For example, 
these properties are mostly all seen in social networks, 
but rarely in technological or biological networks [51]. 
Gronlund and Holme used a model from social network 
analysis to introduce community structure into an Erdos- 
Renyi random network for which clustering, assortativity 
and community structure (using the greedy algorithm) 
were then studied [48] . The evolution of these properties 
was followed as the community structure was introduced. 
Both Q and the clustering coefficient increased from the 
Erdos-Renyi values and showed a strong correlation, fluc- 
tuating at the same times. The networks became assor- 
tative, i.e. there were correlations between degrees of con- 
nected nodes (high degree nodes were more likely to be 
connected to other high degree nodes and low to low). 
The assortativity coefficient fluctuated strongly, but dif- 
ferently to the clustering coefficient and the modularity, 
indicating that, although introducing community struc- 
ture made the networks assortative, the link between the 
two properties is not as strong as that between commu- 
nity structure and clustering. 

Newman studied a model scale-free network with com- 
munity structure [52] that was found to have a high clus- 
tering coefficient. When different sized communities were 
present the model was also assortative by degree [51]. 
The explanation given for this was that nodes with high 
degree are more likely to be in big communities so that 
they can be in the same community as as many of their 
neighbours as possible. Consequently, high degree nodes 



are likely to be in the same communities (the bigger ones) 
and are therefore also likely to be linked to each other. 

The networks studied here have been shown previously 
to be highly clustered [23]. They are also disassortative 
(anticorrelated) by degree, but not significantly differ- 
ently to random networks; this effect can be explained 
by the degree distribution and lack of multiple edges and 
self-connections [43, 47]. The LJ networks have commu- 
nities of differing sizes so they might also be expected 
to be assortative by degree if the hubs are mostly in 
big communities, whereas community sizes in the ran- 
dom networks vary less and as such would be expected 
to be less assortative. As can be seen in plots of the de- 
gree distribution for each community (Figure 8) the hubs 
are split between three of the five communities for LJ13. 
These are not the three largest communities, the highest 
degree node in the second largest community (contain- 
ing 432 nodes) has degree 248 whereas the biggest hub, 
with degree 794, is in a community of only 206 nodes. 
This implies that there is some other property causing 
nodes to group into communities, beyond the degree and 
community sizes. 

The networks are assortative with respect to their po- 
tential energy [24] , i.e. minima with similar potential en- 
ergies are likely to be connected. It seems reasonable 
that nodes grouped in a coininunity would correspond to 
minima with similar energies. This is illustrated in the 
energy distribution shown in Figure 8. There is some sep- 
aration. For example, most of the high energy minima 
are in one community, but the differences are fairly weak. 
In a study of LJ55 a similar overall energy histogram was 
seen and the different peaks could be assigned to different 
classes of structure [53], namely different types and num- 
bers of defect. Previous work on 13-atom Morse clusters 
[28] found that different peaks corresponded to structures 
with different numbers of non-nearest neighbours. It is 
therefore likely that the; different peaks in the energy his- 
togram can also be differentiated by structure. 



B. BLJ13 

A binary Lennard-Jones cluster with 13 atoms, one of 
which is distinguished as being heavy, was also studied. 
The PEL has two funnels [17, 32] leading to different iso- 
mers of the global minimum. Networks for the clusters 
considered in the previous section were found in previous 
work [54] using the eigenvector-following algorithm to lo- 
cate all the stationary points. The network for the binary 
system was obtained from this monoatomic 13-atom net- 
work in the following way. For each minimum in turn, 
one atom was given a mass of two while the other twelve 
had a mass of one. The trace of the inertia tensor was 
calculated as 



f + {Vi - Ucomf + {Zi - Zcomf) (4) 
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FIG. 8: (Colour online) (a) Cumulative degree distribution 
and (b) energy distribution of minima for LJ13. The solid 
(black) line is the overall distribution for the whole network 
and the dashed (coloured) lines those for the five communities 
found using basin hopping. 



where the sum is over all the atoms, rrii is the mass of 
atom i, Xi its x coordinate and Xcom the x-coordinate 
of the centre of mass. This is repeated with a differ- 
ent atom being heavy until all 13 possibilities have been 
investigated. The trace of the inertia tensor is used to 
distinguish different isomers, as the trace of the inertia 
tensor for some of the 13 possible isomers may be equal, 
depending on the symmetry of the cluster. For example, 
the global minimum icosahedron with the heavy atom 
in any of the 12 outside positions are classed as indis- 
tinguishable. There are a total of 17 964 distinguishable 
minima (nodes) compared to 1509 for the monoatomic 
cluster. For each distinguishable version of the transi- 
tion states with one heavy atom, the associated pathway 
was found. The original network tells us which minima 
a transition state connects, and the traces of the inertia 
tensors of the minima at either end of the pathway tells 
us which versions of those minima it connects. This gave 
rise to 294 285 edges in the network. Due to the large 
size of the network, applying the rewiring method to cre- 
ate a random ensemble of networks is not feasible, so the 



original expression for Q has been used. 

The 13-atom cluster with one heavy atom has higher 
modularity than the simple 13-atom cluster, the highest 
Qmax value obtained was 0.4370 using the basin hopping 
algorithm. This reflects the greater heterogeneity of this 
PEL compared to those of the smaller clusters. The two 
permutational isomers of the global minimum are in dif- 
ferent communities, the two largest (consisting of 5507 
and 4732 nodes). Measuring the moment of inertia (the 
geometric mean of the 3 principal moments of inertia) for 
the clusters gives an indication of where the heavy atom 
is; i.e. if it is near the edge of the cluster the moment of 
inertia will be larger and vice versa. The distributions of 
the moments of inertia are shown in Figure 9 for the four 
largest communities. Minima in the largest community, 
containing the version of the global minimum with the 
heavy atom on the outside, have high inertia implying 
that the other minima also have the heavy atom close to 
the outside. The converse is true for the second largest 
community, which contains the permutational isomer of 
the global minimum with the heavy atom in the centre. 
The energy distributions (Figure 9) shows that these two 
communities have very similar energy distributions and 
therefore probably consist of different permutational iso- 
mers of the same geometric structures. There are two 
additional communities which are fairly large (4008 and 
3629 nodes). These also have peaks in their inertia dis- 
tributions at around the same position as the two largest 
communities, implying that the two funnels around the 
two permutational isomers of the global minimum have 
been broken down into smaller communities. One of the 
two smaller communities consists of fairly high energy 
minima and could be a transition region between the two 
funnels. 



C. LJ38 

For LJ38 the community structure is very strong with 
Q'max from the greedy algorithm being 0.8311, which is 
of the order of four times that found for the smaller LJ 
networks and one of the higher modularities found for 
any network. Q' tends to be higher than Q as the lack of 
self-connections means too many edges are generally pre- 
dicted within communities for Q, but Qmax is also greater 
than 0.8. This indicates that the PEL for LJ38 is much 
more heterogeneous than those for the smaller clusters, as 
expected. Q'^ax '^^s also high for the random ensemble 
with a mean of 0.5650 implying that the degree distribu- 
tion can explain some of the community structure. The 
degree distribution of the minima sampled for this clus- 
ter is not scale-free, because of the incompleteness of the 
network. In a scale-free network there are some nodes 
with very high degree (hubs). It is not always likely for 
a hub to be in the same community as all of its neigh- 
bours, for example in LJ14 the largest hub is connected 
to approximately 75 % of the nodes in the network and 
there is no community that large. This could mean that 
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FIG. 9: (Colour online) Distributions of (a) the moments of 
inertia (the geometric mean of the principal moments of in- 
ertia) and (b) energies of minima found for LJ13 with one 
heavy atom, using basin hopping to optimize Q. The solid 
(black) line shows the distribution for all minima. The dashed 
(coloured) lines represent the four largest communities with 
sizes greater than 3000, the other communities all contained 
less than 20 nodes. Minima with larger moments of inertia 
have the heavy atom closer to the outside of the cluster. 



community structure can be stronger in networks with 
degree distributions without hubs, such as that for LJ38, 
explaining the high values of the modularity seen for the 
random ensemble. However, the value of Q'max fo'" the 
PEL network is over 17 standard deviations higher than 
that for the random networks and so is significant. The 
dendrogram is shown in Figure 10, only the last 100 of 
6000 steps are shown. 

The community structure found can be compared to 
the multiple-funnel structure of the PEL. Some of the 
minima have previously been assigned to either an fee 
or icosahedral funnel using a master equation approach. 
At the point of maximum Q' in the community structure 
found from the greedy algorithm, the minima from the fee 
funnel are in the same community and there are no icosa- 
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FIG. 10: The final 100 steps in the dendrogram for LJ38 using 
the greedy algorithm. If a community only consists of single 
nodes, i.e. is not composed of smaller communities, then for 
clarity only the first and last nodes to join that community 
are shown. The variation of Q' as the algorithm progresses is 
not shown, it has a single peak of value (3^^1=0.8311 at level 
5956 (indicated by a dashed line on the dendrogram). 



hedral structures in that community. This is also true for 
the partition with the highest value of Q' ((5'=0.8314), 
obtained from basin hopping. The best split is at level 
5956 ((5^„j.=0.8311) and the fee and icosahedral min- 
ima are in different communities until level 5990 where 
(5'=0.8099 and has begun to fall off. The betweenness 
algorithm is too slow to run to completion, but after re- 
moval of only 13 edges the network was split into two 
communities with all fee minima in one and all icosahe- 
dral minima in the other. The community algorithms 
are therefore giving an insight into the topology of the 
PEL whilst only using data on the connectivities of the 
minima, whereas the previous methods also use the en- 
ergies of minima and the barrier heights between them. 
This should be contrasted with the community structure 
found for the smaller clusters, which is very weak and 
was not found at all using the betweenness algorithm. 
The PELs for the small clusters have a single funnel to- 
pography. The 13-atom Morse cluster, which is similar to 
LJ13 when appropriate parameters are used, has only one 
monotonic sequence basin, meaning that a path from any 
minimum which decreases the energy at each step leads 
to the global minimum [29] . This information from topo- 
graphical analysis is consistent with the picture obtained 
from the topological approach in this work. 



D. 2D hexagonal lattice 

Since the energy landscapes of small clusters have sin- 
gle funnel topographies, it is curious that there should 
be any community structure in the corresponding net- 
works. It is possible that this community structure is 
related to the spatial nature of the landscapes. Nodes 
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in the same region of space are more likely to be inter- 
connected (leading to high clustering). If these nodes are 
grouped together into communities, there is likely to be 
many edges within those communities and few between 
different communities, as edges connecting different com- 
munities would only come from nodes adjacent to the 
boundaries. To investigate the effect of spatial organiza- 
tion and clustering on community structure a hexagonal 
lattice was studied. 

The greedy algorithm can be applied to a large two- 
dimensional hexagonal lattice, where all nodes have the 
same degree (fc=6) and the clustering coefficient is 0.4. 
In the first step any edge is equally likely to be added, 

each changing Q to Q -I- AQ, where AQ = If 
P/2M < 1, i.e. M > 18 then AQ is positive. In the sec- 
ond step one of the two nodes that forms a triangle with 
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that edge will be added, with AQ = 2 x 



(l-fc^/2M) 
M 



(also 



positive if AI > 18). In the third step, a fourth node will 
be added to the community. This can only increase the 
number of edges within the community by two, whereas 
the increase in the number of expected edges will be the 
number of new pairs of nodes (the size of the commu- 
nity, Uc = 3) multiplied by the probability of an edge 
between each pair, AQ = 2-"cfc^/2M ^ ^j^^ following 
steps, the number of edges added to the community can 
be increased by up to 6 (the degree of each node) depend- 
ing on the shape of the growing community (this will de- 
pend on the random choices). At each step the expected 
number of edges (to be subtracted) is Uck'^ /2M. At some 
point, AQ will be smaller than that for starting a new 

community (AQ = ^ ^^^^ community will 

be started. 

Q and Q' were optimized for a hexagonal lattice with 
periodic boundary conditions such that all nodes have 
degree 6. The following results correspond to the largest 
system studied, a lattice with 2500 nodes. These were 
also compared to a random ensemble. Qmax and Q'^^ax 
from greedy optimization are 0.7674 and 0.7231 respec- 
tively, whereas the corresponding values for the random 
networks are 0.3860 and 0.3517, differences of 191 and 
45 standard deviations. The highest value of Q' (0.8480) 
corresponds to the partition shown in Figure 11, which 
was found by basin hopping from a partition into roughly 
equally sized hexagons (Q'=0.8389 for 18 communities of 
average size 139). Interestingly, the partition obtained 
from greedy optimization contained much larger commu- 
nities, with an average size over 300. Nodes in the lattice 
have constant degree, so this could possibly explain why 
the modularity is higher than for scale- free landscape net- 
works (as for LJas), for both the lattice and the random- 
ized versions. However, Apollonian networks, which are 
spatial scale-free networks, also have a high modularity. 
The modularity of a scale-free network embedded on a 
lattice could be studied [15] to provide a further compar- 
ison. 

The modularity of the lattice is still much higher than 
that for the random ensemble so the community structure 



"o*o*o* 

'o''o''o"o''o''o"o"'o''o''o"o o o o o''a"'o''o''o''o'' 

■oVoVoVoVoVoSVoSSSV." 
'oVoSSVoV.VoV.VoSW 

looaaooaeeoaaeoaee 
I o o o o o o o o o d o o 

oSVoSSVo 

v.^. oaeoeeeoo 



•oWoW.y 

.WoV 



'oVo° 
) o o I 



oooooooao 
eaoaeeooea 

■oVoSVoVoSVo 

1 o.o_o_o O_0_0_O O_0_( 



o o o o o 



ooaoooaoooaaoooc 
iQaaoooaoooooooo 
oooeeQaaoeooooa< 
..-.jQoeeeooeoaaeeea 
eooaeeeeeeoaeoaeeee< 
^ n Q O O O 0_0 OOaoOOooOOa 
o a e a o o o o ooooooooc 



ooooooaaeooeao 
_ ...eoooeeoeeeoaeet 



loooaoooaoeo 



^QOOOOOOOQDOOOOOQOOI 

aaaaeQoaoooaeoQeeeo 
>oaaeeeaeeeaeeeoeeoi 



eooeoQee 
'oVoVoVoV. 

WoVoVoV. 

'q o o Q^o^a 0°c 

looaoeaee 

eoooeaoeeoaaeoof 



loaooaoooa 

.VoVoVoVoV 



) o o o o o o o o o o o < 



looeeooeoooeeooeeooee 



Qoooooooeo 
looeeooeee 

ooooaaaoof 

OOOOOOOOOI 

oooooooeoe 



FIG. 11; (Colour online) The partition with the highest mod- 
ularity (Q'=0.8480) seen for a two-dimensional hexagonal lat- 
tice with 2500 nodes using periodic boundary conditions. This 
partition was found using basin hopping to optimize Q' from 
an initial partition into roughly equally sized hexagons. 



is significant. Nodes close together on a hexagonal lattice 
have many edges between them due to the highly clus- 
tered nature of the network. They also have few edges to 
the rest of the network, as they only occur at the edges 
of the community. Similar high values of the modular- 
ity have been analytically predicted for low-dimensional 
regular lattices [49] . The partition of the lattice into com- 
munities is highly degenerate, for example the partitions 
could be translated in any direction, giving a different 
community split with the same high value of the modu- 
larity. This raises questions about the interpretation of 
Q and Qrand in terms of 'unique' communities. Other 
methods for finding community structure have been pro- 
posed which give many possible partitions due to random 
choices in the algorithms [35, 46]. These partitions could 
then be compared to determine whether they are based 
on a strong community structure or are very different to 
each other, as would be the case for the hexagonal lattice. 

An alternative way to probe this issue would be to 
look at the thermodynamics associated with partitioning 
of the network, where — Q or — Q' plays the role of energy, 
as in the Monte Carlo methods used here. If there was 
a unique, well-defined community structure, one would 
expect a transition, somewhat akin to crystallization, to 
a low-temperature state with a low density of partitions, 
all of which would have the same basic structure. By 
contrast, a more 'glass-like' behaviour would be expected 
for networks with ill-defined community structure. That 
is there would only be a gradual increase in Q as the 
temperature is decreased without any sharp transitions, 
and there would be many high-modularity partitions with 
significant structural differences. 



V. CONCLUSION 

The potential energy landscapes of small Lennard- 
Jones clusters have been studied in terms of networks 
describing which metastable states of the clusters are 
connected. A number of algorithms have been used to 
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uncover the community structure of the PEL networks. 
Ah are based solely on the topology of the network, 
i.e. barrier heights and energies of the minima are not 
taken into account. The first algorithm used, based on 
edge betweenness, found no community structure in the 
networks of the small clusters, consistent with a single 
funnel shape of the PEL. However, this algorithm has 
been unsuccessful in previous studies when the network 
is densely connected, and the PEL networks studied here 
have a higher average degree than many of the other 
networks previously studied. The other algorithms put 
nodes into commimities such that there are more edges 
within communities than expected for a random net- 
work, thus optimizing the modularity. A modification 
was made to the way in which the community structure 
is quantified. Q' is essentially the same as Q, measur- 
ing the fraction of edges within communities compared 
to the predicted fraction for a random network with the 
same degree distribution. However, here the random net- 
work does not contain multiple edges or self-connections 
making it a better comparison. 

Optimizing the modularity did find weak community 
structure in the networks. This implies that within the 
single funnel of the PEL minima form groups with more 
transition states between them than to minima in other 
groups. The community structure is likely to be con- 
nected to the spatial nature of the networks and therefore 



related to the clustering. To investigate the effect of spa- 
tial ordering of a network on its community structure a 
hexagonal lattice was studied. Strong community struc- 
ture was found for the lattice, implying that the weak 
community structure found in the landscape networks 
could be due to spatial organization. 

Two LJ clusters with more complicated PELs where 

community structure was expected were also studied. La- 
belling one atom in LJ13 distinguishes permutational iso- 
mers that fit broadly into two classes, those with the 
heavy atom on the outside and those with the heavy 
atom on the inside. Two communities corresponding to 
these classes of structure were found by optimizing Q. 
The PEL of LJ38 is known to have two funnels, one con- 
sisting of fee structures containing the truncated octahe- 
dron global minimum and the other consisting of icosahe- 
dral structures. Both topological algorithms found much 
stronger community structure for this cluster than for 
the smaller clusters, implying a much more heterogeneous 
PEL. The community structure found was also compared 
to the known funnel structure of the PEL. In both algo- 
rithms fee and icosahedral structures were in different 
communities. The topology of the PEL therefore con- 
tains information about its heterogeneity and can be used 
to provide a similar picture to that gained from a topo- 
graphical analysis. 
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